Datos

En esta práctica exploraremos una aplicación en astronomía. Los datos los hemos tomado de aquí y buscaremos el mismo objetivo: predecir el diámetro de un asteroide en términos de ciertas características medidas por instrumentos satelitales y otros sensores.

La base de datos proviene del JPL (Jet Propulsion Laboratory) en Caltech y pueden buscar los datos crudos aquí. Nosotros usaremos una base de datos previamente procesada que contiene predictores para esta tarea.

Los datos tienen severos problemas de valores atípicos en la escala original de diametro (medido en metros). Esto es consecuencia de la regla de potencias descrita en la imagen de arriba. Para realizar la predicción consideraremos utilizar los datos en escala logaritmica.

Asimismo, algunos de los predictores podrían necesitar un tratamiento similar, pues algunos también incluyen valores atípicos y/o son atributos que sólo presentan valores positivos.

Objetivo

El objetivo es construir modelos predictivos, compararlos y escoger el mejor de ellos. Hemos desarrollado el código sobre una máquina con las siguientes características de sistema operativo, CPU (procesador y núcleos) y RAM respectivamente:

## [1] "macOS Monterey 12.2"

## [1] "Apple M1"

## [1] 8

## 8.6 GB

Por último, las predicciones se harán para un conjunto de prueba para el cual no tienen acceso a la variable objetivo. La métrica de evaluación será el error cuadrático medio en escala logarítmica, es decir,

\[ \mathscr{L}_V(y, \hat y ) = \frac1m \sum \left(\log y_i - \log \hat y_i \right)^2\,.\]

Por lo tanto, se adjuntará al proyecto un archivo *.csv con las predicciones en escala logarítmica: \(\log(\text{diametro})\).

Exploración y limpieza

En primer lugar, buscamos valores faltantes.

na_values <- c()
for (i in 1:ncol(train)) {
  na_values[i] <- sum(is.na(train[, i]))
}
na_values

##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Al no haber datos faltantes, no necesitamos imputar por ningún método.

Como segundo paso, observamos que la mayoría de los atributos son numéricos, excepto por el atributo class.

count(train, class)

## # A tibble: 9 × 2
##   class     n
##   <chr> <int>
## 1 AMO      30
## 2 APO      44
## 3 ATE      12
## 4 CEN       3
## 5 IMB      66
## 6 MBA   12497
## 7 MCA      33
## 8 OMB     732
## 9 TJN     184

Es fácil ver que una gran mayoría de los asteroides observados son de la clase MBA tal que usar este atributo para nuestros modelos de predicción se vuelve inútil. Por lo tanto, redefinimos nuestros conjuntos de prueba y entrenamiento a escala logarítmica como sigue:

test_1 <- log(subset(test, select = -c(class)))
train_1 <- log(subset(train, select = -c(class)))

El atributo condition_code se vuelve inútil después de convertirlo a la forma log(condition_code) pues toma valores desde \(-\infty\). Nótese, por el mismo efecto, su correlación con el diámetro: NaN.

train_cor <- as.data.frame(cor(train_1))
diameter_cor <- subset(train_cor, select = c(diameter))

diameter_cor %>% arrange(abs(diameter))

##                diameter
## w                0.0062
## om               0.0120
## ma               0.0210
## i                0.0962
## e               -0.1650
## albedo          -0.2504
## data_arc         0.3109
## n_obs_used       0.4493
## ad               0.4809
## moid             0.5217
## q                0.5419
## a                0.5680
## n               -0.5680
## per              0.5680
## per_y            0.5680
## H               -0.8430
## diameter         1.0000
## condition_code      NaN

Redefinimos nuestros conjuntos como:

test_set <- subset(test_1, select = -c(condition_code))
train_set <- subset(train_1, select = -c(condition_code))

Training

Normalizamos todos los datos a través de nuestra receta para mejorar las predicciones de los modelos.

receta <- recipe(diameter ~ ., train_set) %>% 
  step_normalize(all_predictors()) %>% 
  prep()

Todos los modelos siguen la siguiente estructura:

modelo_n, definimos el modelo a utilizar libre de parámetros;
vc, dividimos nuestro conjunto de entrenamiento en bloques (por vfold_cv() o bootstrap());
flujo_n, creamos el workflow con el modelo_n y nuestro nuestra receta;
grid_n, definimos nuestro grid de búsqueda de (hiper-)parámetros (por grid_latin_hypercube(), grid_random() o de forma manual como secuencia);
metricas_n, ejecutamos el proceso de ajuste de (hiper-)parámetros;
mejor_n, seleccionamos las mejores métricas son select_best();
modelo_final_n, redefinimos el modelo con los parámetros óptimos;
flujo_final_n, creamos el workflow final con el modelo_final_n y nuestro nuestra receta;
ajuste_final_n, a través de nuestro flujo_final_n ajustamos nuestro modelo a nuestro conjunto de entrenamiento;
pred_n, sacamos nuestras predicciones del mismo conjunto de entrenamiento con nuestro modelo ajustado;
p_n, graficamos un histograma para comparar nuestras predicciones con las etiquetas reales;
q_n, hicimos un scatter plot de etiquetas reales contra nuestras predicciones;
score[n], determinamos el coeficiente de determinación de cada una de las predicciones para escoger el mejor modelo a través de un score.

Regresión lineal regularizada

Ridge

# Modelo
modelo_1 <- linear_reg(mixture = 0, penalty = 0) %>%
  set_engine("glmnet") %>%
  set_mode("regression") %>%
  fit(diameter ~ ., juice(receta))

coefs_1 <- tidy(modelo_1$fit) %>%
  filter(term != "(Intercept)")

Realizamos este plot para cerciorarnos que nuestra penalización \(\lambda\) converja.

Validación cruzada y ajuste de parámetros

Por lo observado en la gráfica anterior, podemos optar con seguridad por un rango de penalización entre 1e-01 y 1e+04.

modelo_1_regularizado <-
  linear_reg(mixture = 0, penalty = tune()) %>%
  set_engine("glmnet")

flujo_1 <- workflow() %>%
  add_model(modelo_1_regularizado) %>%
  add_recipe(receta)

bf_set_1 <-
  parameters(penalty(range = c(-1, 4), trans = log10_trans()))

grid_1 <- grid_regular(bf_set_1, levels = 50)

vc <- vfold_cv(train_set)

tic()
metricas_1 <- tune_grid(
  flujo_1,
  resamples = vc,
  grid = grid_1,
  metrics = metric_set(rmse)
)
toc()

## 4.502 sec elapsed

mejor_1 <- metricas_1 %>% select_best()

Visualización de desempeño

Nótese la baja variabilidad en el error de validación tal que, después de pruebas con intervalos más grandes al definido en el subapartado anterior, \(\lambda \to 0\). Esto se debe a que anteriormente, hicimos una transformación logarítmica a nuestro conjunto de datos orignal.

Predicciones

modelo_final_1 <- finalize_model(modelo_1_regularizado, mejor_1)

flujo_final_1 <-
  workflow() %>% add_model(modelo_final_1) %>% add_recipe(receta)

tic()
ajuste_final_1 <- flujo_final_1 %>% fit(train_set)
toc()

## 0.186 sec elapsed

pred_1 <- predict(ajuste_final_1, new_data = train_set)

score[1] <-
  R2_Score(y_pred = pred_1$.pred, y_true = train_set$diameter)

La puntuación \(R^2\) obtenida con esta regresión es 0.94.

Lasso

# Modelo
modelo_2 <- linear_reg(mixture = 1, penalty = 0) %>%
  set_engine("glmnet") %>%
  set_mode("regression") %>%
  fit(diameter ~ ., juice(receta))

coefs_2 <- tidy(modelo_2$fit) %>%
  filter(term != "(Intercept)")

Realizamos este plot para cerciorarnos que nuestra penalización \(\lambda\) converja.

Validación cruzada y ajuste de parámetros

Por lo observado en la gráfica anterior, podemos optar con seguridad por un rango de penalización entre 1e-02 y 1e+01. Sin embargo, la convergencia no parece ser la deseada.

modelo_2_regularizado <-
  linear_reg(mixture = 1, penalty = tune()) %>%
  set_engine("glmnet")

flujo_2 <- workflow() %>%
  add_model(modelo_2_regularizado) %>%
  add_recipe(receta)

bf_set_2 <-
  parameters(penalty(range = c(-2, 1), trans = log10_trans()))

grid_2 <- grid_regular(bf_set_2, levels = 50)

vc <- vfold_cv(train_set)

tic()
metricas_2 <- tune_grid(
  flujo_2,
  resamples = vc,
  grid = grid_2,
  metrics = metric_set(rmse)
)
toc()

## 4.786 sec elapsed

mejor_2 <- metricas_2 %>% select_best()

Visualización de desempeño

En este caso, se puede observar el mismo efecto del apartado anterior, pues nuestro error de validación \(\lambda \to 0\).

Predicciones

modelo_final_2 <- finalize_model(modelo_2_regularizado, mejor_2)

flujo_final_2 <-
  workflow() %>% add_model(modelo_final_2) %>% add_recipe(receta)

tic()
ajuste_final_2 <- flujo_final_2 %>% fit(train_set)
toc()

## 0.202 sec elapsed

pred_2 <- predict(ajuste_final_2, new_data = train_set)

score[2] <-
  R2_Score(y_pred = pred_1$.pred, y_true = train_set$diameter)

La puntuación \(R^2\) obtenida con esta regresión es 0.94.

Elastic net

# Modelo
modelo_3 <-
  linear_reg(mixture = tune(), penalty = tune()) %>%
  set_engine("glmnet")

Validación cruzada y ajuste de parámetros

flujo_3 <- workflow() %>%
  add_model(modelo_3) %>%
  add_recipe(receta)

bf_set_3 <-
  parameters(penalty(range = c(-2, 2), trans = log10_trans()), mixture(range = c(0, 1)))

grid_3 <- grid_regular(bf_set_3, levels = 50)

vc <- vfold_cv(train_set)

tic()
metricas_3 <- tune_grid(
  flujo_3,
  resamples = vc,
  grid = grid_3,
  metrics = metric_set(rmse),
  control = control_grid(allow_par = TRUE) # Aquí habilitamos el cómputo en paralelo.
)
toc()

## 185.223 sec elapsed

mejor_3 <- metricas_3 %>% select_best()

Visualización de desempeño

En este caso, después de algunas pruebas, nos decidimos por difinir el intervalo de penalización de 1e-02 a 1e+02.

Predicciones

modelo_final_3 <- finalize_model(modelo_3, mejor_3)

flujo_final_3 <-
  workflow() %>% add_model(modelo_final_3) %>% add_recipe(receta)

tic()
ajuste_final_3 <- flujo_final_3 %>% fit(train_set)
toc()

## 0.189 sec elapsed

pred_3 <- predict(ajuste_final_3, new_data = train_set)

score[3] <-
  R2_Score(y_pred = pred_3$.pred, y_true = train_set$diameter)

La puntuación \(R^2\) obtenida con esta regresión es 0.96.

Árboles de decisión

# Modelo
modelo_4 <-
  decision_tree(cost_complexity = tune(), tree_depth = tune()) %>%
  set_engine("rpart") %>%
  set_mode("regression")

Validación cruzada y ajuste de parámetros

grid_4 <- grid_regular(cost_complexity(),
                          tree_depth(),
                          levels = 5)

vc <- vfold_cv(train_set)

tic()
metricas_4 <- tune_grid(
  modelo_4,
  diameter ~ .,
  resamples = vc,
  grid = grid_4,
  metrics = metric_set(rmse),
  control = control_grid(allow_par = TRUE) # Aquí habilitamos el cómputo en paralelo.
)
toc()

## 68.288 sec elapsed

mejor_4 <- metricas_4 %>% select_best()

Predicciones

modelo_final_4 <- finalize_model(modelo_4, mejor_4)

flujo_final_4 <-
  workflow() %>% add_model(modelo_final_4) %>% add_recipe(receta)

tic()
ajuste_final_4 <- flujo_final_4 %>% fit(train_set)
toc()

## 0.684 sec elapsed

pred_4 <- predict(ajuste_final_4, new_data = train_set)

score[4] <-
  R2_Score(y_pred = pred_4$.pred, y_true = train_set$diameter)

La puntuación \(R^2\) obtenida con esta regresión es 0.98.

Bosques aleatorios

# Modelo
modelo_5 <-
  rand_forest(mtry = tune(),
              trees = tune(),
              min_n = tune()) %>%
  set_engine("ranger") %>%
  set_mode("regression")

Validación cruzada y ajuste de parámetros

vc <- bootstraps(train_set, times = 5)

flujo_5 <-
  workflow() %>%
  add_recipe(receta) %>%
  add_model(modelo_5)

grid_5 <- parameters(
   finalize(mtry(), select(train_set, -diameter)), # Utilizamos la función finalize()
                                                   # para que tardara menos y evitar un warning.
   trees(),
   min_n()) %>% 
   grid_random(5) # A diferencia de los *grid* de búsqueda de los demás modelos,
                  # en esta ocasión lo definimos de forma aleatoria.

tic()
metricas_5 <- tune_race_anova(
  flujo_5,
  resamples = vc,
  grid = grid_5,
  metrics = metric_set(rmse),
  control = control_race(allow_par = TRUE) # Aquí habilitamos el cómputo en paralelo.
)
toc()

## 361.568 sec elapsed

mejor_5 <- metricas_5 %>% select_best()

Predicciones

modelo_final_5 <- finalize_model(modelo_5, mejor_5)

flujo_final_5 <-
  workflow() %>% add_model(modelo_final_5) %>% add_recipe(receta)

tic()
ajuste_final_5 <- flujo_final_5 %>% fit(train_set)
toc()

## 100.299 sec elapsed

pred_5 <- predict(ajuste_final_5, new_data = train_set)

score[5] <-
  R2_Score(y_pred = pred_5$.pred, y_true = train_set$diameter)

La puntuación \(R^2\) obtenida con esta regresión es 1.

Máquinas de soporte vectorial

Elegimos un modelo de máquinas de soporte vectorial con función de base radial.

# Modelo
modelo_6 <-
  svm_rbf(cost = tune(),
          rbf_sigma = tune()) %>%
  set_engine("kernlab") %>%
  set_mode("regression")

Validación cruzada y ajuste de parámetros

vc <- bootstraps(train_set, times = 5)

flujo_6 <-
  workflow() %>%
  add_recipe(receta) %>%
  add_model(modelo_6)

# Definimos nuestro grid de búsqueda como un hipercubo de tamaño n = 5.
grid_6 <- modelo_6 %>% parameters() %>% grid_latin_hypercube(size = 5)

tic()
metricas_6 <- tune_race_anova(
  flujo_6,
  resamples = vc,
  grid = grid_6,
  metrics = metric_set(rmse),
  control = control_race(allow_par = TRUE) # Aquí habilitamos el cómputo en paralelo.
)
toc()

## 159.674 sec elapsed

mejor_6 <- metricas_6 %>% select_best()

Predicciones

modelo_final_6 <- finalize_model(modelo_6, mejor_6)

flujo_final_6 <-
  workflow() %>% add_model(modelo_final_6) %>% add_recipe(receta)

tic()
ajuste_final_6 <- flujo_final_6 %>% fit(train_set)
toc()

## 20.554 sec elapsed

pred_6 <- predict(ajuste_final_6, new_data = train_set)

score[6] <-
  R2_Score(y_pred = pred_6$.pred, y_true = train_set$diameter)

La puntuación \(R^2\) obtenida con esta regresión es 0.85.

XGBoost

# Modelo
modelo_7 <- boost_tree(
  trees = 1000,
  tree_depth = tune(),
  min_n = tune(),
  loss_reduction = tune(),
  sample_size = tune(),
  mtry = tune(),
  learn_rate = tune(),
) %>%
  set_engine("xgboost") %>%
  set_mode("regression")

Validación cruzada y ajuste de parámetros

vc <- vfold_cv(train_set, strata = diameter)

grid_7 <- grid_latin_hypercube(  # Definimos nuestro grid de búsqueda como un
                                 # hipercubo de tamaño 10.
  tree_depth(),
  min_n(),
  loss_reduction(),
  sample_size = sample_prop(),
  finalize(mtry(), train_set),
  learn_rate(),
  size = 10
)

flujo_7 <- workflow() %>%
  add_recipe(receta) %>%
  add_model(modelo_7)

tic()
metricas_7 <- tune_race_anova(
  flujo_7,
  resamples = vc,
  grid = grid_7,
  metrics = metric_set(rmse),
  control = control_race(allow_par = TRUE, parallel_over = "resamples") 
) # Habilitamos el cómputo en paralelo solo en el remuestreo
  # para reducir el tiempo de ejecución.
toc()

## 209.883 sec elapsed

mejor_7 <- metricas_7 %>% select_best()

Predicciones

modelo_final_7 <- finalize_model(modelo_7, mejor_7)

flujo_final_7 <-
  workflow() %>% add_model(modelo_final_7) %>% add_recipe(receta)

tic()
ajuste_final_7 <- flujo_final_7 %>% fit(train_set)
toc()

## 24.299 sec elapsed

pred_7 <- predict(ajuste_final_7, new_data = train_set)

score[7] <-
  R2_Score(y_pred = pred_7$.pred, y_true = train_set$diameter)

La puntuación \(R^2\) obtenida con esta regresión es 0.98.

Importancia de las variables

## Warning: `pull_workflow_fit()` was deprecated in workflows 0.2.3.
## Please use `extract_fit_parsnip()` instead.

Modelo final

Conclusiones

A pesar de que no fijamos una semilla para la aleatoriedad, el modelo de Bosques Aleatorios obtuvo siempre el mejor desempeño y las mejores predicciones. Esto nos conduce a pensar que no existe multicolinealidad entre los atributos. Por lo tanto, nuestras predicciones en escala logarítmica son las siguientes:

modelo_final <- ajuste_final_5
write_csv(as.data.frame(predict(modelo_final, new_data = test_set)$.pred),
          'predicciones_04.csv')

Proyecto final

Aprendizaje estadístico

Carlos Lezama, Alejandro Chávez, Jorge Rizo, Manuel Garduño

mayo 2021

Datos

Objetivo

Exploración y limpieza

Training

Regresión lineal regularizada

Ridge

Validación cruzada y ajuste de parámetros

Visualización de desempeño

Predicciones

Lasso

Validación cruzada y ajuste de parámetros

Visualización de desempeño

Predicciones

Elastic net

Validación cruzada y ajuste de parámetros

Visualización de desempeño

Predicciones

Árboles de decisión

Validación cruzada y ajuste de parámetros

Predicciones

Bosques aleatorios

Validación cruzada y ajuste de parámetros

Predicciones

Máquinas de soporte vectorial

Validación cruzada y ajuste de parámetros

Predicciones

XGBoost

Validación cruzada y ajuste de parámetros

Predicciones

Importancia de las variables

Modelo final

Conclusiones